# 文件创建时间: 2025-09-24
# 文件功能描述: 两相流物理模型类，包括流体物性计算、壁面摩擦、界面摩擦和湍流混合等物理效应

from grid import SubChannelType
import numpy as np
from CoolProp.CoolProp import PropsSI
from enum import Enum
from typing import Tuple, Dict

class PhaseType(Enum):
    LIQUID = "liquid"
    VAPOR = "vapor"

class PhysicsModel:
    """两相流物理模型类"""
    def __init__(self, fluid_properties: Dict[str, float]):
        # 封闭关系常数
        self.BETA = 0.065      # 湍流交混系数
        self.FTM = 0.8         # 端流动量因子
        self.KIJ = 0.5         # 间隙横向流阻
        self.C_VM = 20.0       # 虚拟质量系数
        self.C_i = 1.0         # 界面阻力系数（【修复2025-11-04】从0.3增加到1.0，加强相间耦合）
        self.lateral_flow_resistance = 0.7  # 横向流动阻力系数（默认值，会被配置覆盖）
        
        # 子通道类型系数
        self.channel_coefficients = {
            SubChannelType.CORNER: {
                'wall_friction': 1.2,    # 角通道壁面摩擦增大
                'lateral_flow': 0.8,     # 角通道横向流动减小
                'turbulent_mixing': 0.7  # 角通道湍流混合减小
            },
            SubChannelType.EDGE: {
                'wall_friction': 1.1,    # 边通道壁面摩擦略增
                'lateral_flow': 1.0,     # 边通道横向流动正常
                'turbulent_mixing': 1.0  # 边通道湍流混合正常
            },
            SubChannelType.CENTER: {
                'wall_friction': 1.0,    # 中心通道壁面摩擦正常
                'lateral_flow': 1.2,     # 中心通道横向流动增大
                'turbulent_mixing': 1.1  # 中心通道湍流混合增大
            }
        }
        
        # 重力加速度
        self.g = 9.81          # m/s^2
        
        # TODO: 任务 3.4 - 参数化输入
        # 从配置中获取流体物性
        self.T_ref = fluid_properties.get('temperature', 298.15)
        self.liquid_density = fluid_properties.get('liquid_density', 997.0)
        self.liquid_viscosity = fluid_properties.get('liquid_viscosity', 8.9e-4)
        self.vapor_density = fluid_properties.get('vapor_density', 1.225)
        self.vapor_viscosity = fluid_properties.get('vapor_viscosity', 1.81e-5)

        
        # 物性缓存
        self._properties_cache = {}
        self._pressure_range = (0.5e5, 2e5)  # 压力范围：0.5-2 bar
        
        # 初始化流型识别器
        from flow_regimes import FlowRegimeIdentifier
        self.flow_regime_identifier = FlowRegimeIdentifier()
        
    def _get_cache_key(self, P, phase):
        """生成缓存键，按kPa取整"""
        P_kpa = round(P / 1000) * 1000
        return (P_kpa, phase)
    
    def calculate_properties(self, P, phase):
        """计算流体物性，支持标量或数组输入
        Args:
            P: 压力 [Pa] (可以是标量或NumPy数组)
            phase: 相态 (液相或气相)
        Returns:
            properties: 包含密度、粘度的字典，值为标量或NumPy数组
        """
        is_array = isinstance(P, np.ndarray)

        if phase == PhaseType.LIQUID:
            density = self.liquid_density
            viscosity = self.liquid_viscosity
        else:  # VAPOR (空气)
            density = self.vapor_density
            viscosity = self.vapor_viscosity

        if is_array:
            return {
                'density': np.full_like(P, density),
                'viscosity': np.full_like(P, viscosity)
            }
        else:
            return {
                'density': density,
                'viscosity': viscosity
            }
    
    def calculate_wall_friction(self, phase_type, density, velocity, d_h, viscosity, alpha, channel_type):
        """计算壁面摩擦力
        
        根据两相流摩擦模型计算壁面摩擦力，参考COBRA软件的摩擦力模型
        
        Args:
            phase_type: 相态（气相或液相）
            density: 密度
            velocity: 速度
            d_h: 水力直径
            viscosity: 动力粘度
            alpha: 空泡份额
            channel_type: 子通道类型
        Returns:
            tau_w: 壁面摩擦力
        """
        # 保护参数，防止无效值导致计算不稳定
        # 防止无效值的有效性检查
        if not (np.isfinite(density) and np.isfinite(velocity) and np.isfinite(d_h) and 
                np.isfinite(viscosity) and np.isfinite(alpha)):
            return 0.0
            
        velocity_abs = abs(velocity)
        alpha_safe = np.clip(alpha, 1e-6, 1.0-1e-6)
        density_safe = max(density, 1e-6)
        viscosity_safe = max(viscosity, 1e-10)
        d_h_safe = max(d_h, 1e-6)
        
        # 计算雷诺数 - 不将空泡率纳入雷诺数计算，这样更符合物理含义
        Re = density_safe * velocity_abs * d_h_safe / viscosity_safe
        
        # 计算基本摩擦系数 - 改进摩擦系数计算
        if Re < 2000:
            # 层流区域
            f = 64.0 / max(Re, 1.0)  
        elif Re < 4000:
            # 过渡区域 - 平滑过渡
            f_lam = 64.0 / max(Re, 1.0)
            f_turb = 0.316 / (max(Re, 1.0)**0.25)
            ratio = (Re - 2000) / 2000.0  # 过渡比例
            f = f_lam * (1 - ratio) + f_turb * ratio
        else:
            # 湍流区域 - 使用Blasius公式
            f = 0.316 / (max(Re, 1.0)**0.25)
        
        # 限制摩擦系数最大值，提高数值稳定性
        f = min(f, 0.1)  
        
        # 子通道类型修正 - 保留原有修正
        f *= self.channel_coefficients[channel_type]['wall_friction']
        
        # 壁面接触面积分数 - 基于流型的分段函数，参考COBRA软件
        if phase_type == PhaseType.VAPOR:
            # 气相壁面接触分数
            if alpha_safe < 0.3:
                # 气泡流区域 - 接触面积小
                contact_fraction = 0.05 * alpha_safe**0.5
            elif alpha_safe < 0.8:
                # 过渡区域 - 接触面积适中
                contact_fraction = 0.2 * alpha_safe**0.3
            else:
                # 环状流区域 - 接触面积大
                contact_fraction = 0.5 * alpha_safe**0.2
        else:  # 液相
            # 液相壁面接触分数
            if alpha_safe < 0.3:
                # 单相液体区域 - 接触面积最大
                contact_fraction = 0.95 - 0.2 * alpha_safe
            elif alpha_safe < 0.8:
                # 过渡区域
                contact_fraction = 0.8 * (1-alpha_safe)**0.3
            else:
                # 液膜区域 - 接触面积小
                contact_fraction = 0.6 * (1-alpha_safe)**0.8
        
        # 确保接触面积系数合理
        contact_fraction = np.clip(contact_fraction, 0.01, 0.99)
        
        # 计算两相流摩擦力乘数 - 参考COBRA的Armand相关式
        if phase_type == PhaseType.VAPOR:
            if alpha_safe <= 0.6:
                phi_2 = 1.0  # 气相在低空泡率时基本不受影响
            else:
                # 气相在高空泡率时摩擦系数增大
                phi_2 = 1.0 + 0.5 * (alpha_safe - 0.6) / 0.4
        else:  # 液相
            if alpha_safe <= 0.6:
                # Armand相关式简化版
                phi_2 = (1.0 - alpha_safe)**(-1.42)
            else:
                # 高空泡率时液相摩擦减小
                phi_2 = 0.5 * (1.0 - alpha_safe)**(-1.0)
        
        # 限制摩擦乘数避免极端值
        phi_2 = min(phi_2, 10.0)
        
        # 计算壁面摩擦力 - 优化：降低系数以提高数值稳定性
        tau_w = -0.03 * f * density_safe * velocity_abs * velocity * contact_fraction * phi_2
        
        # 低速区线性化处理
        if velocity_abs < 0.05:
            tau_w = -0.03 * f * density_safe * 0.05 * velocity * contact_fraction * phi_2
        
        # 摩擦力限幅，防止极端值（优化：从10.0降低到5.0）
        max_friction = 5.0
        if abs(tau_w) > max_friction:
            tau_w = -max_friction * np.sign(velocity)
        
        # 最终有效性检查
        if not np.isfinite(tau_w):
            return 0.0
        
        return tau_w
        
    def calculate_interfacial_friction(self, alpha, rho_g, rho_l, v_g, v_l, d_h, channel_type, flow_regime=None, mu_g=1.8e-5, mu_l=1.0e-3, temperature=298.15, pressure=1.0e5):
        """计算界面剪切力
        
        Args:
            alpha: 空泡份额
            rho_g: 气相密度
            rho_l: 液相密度
            v_g: 气相速度
            v_l: 液相速度
            d_h: 水力直径
            channel_type: 子通道类型
            flow_regime: 流型字符串(可选)，如果提供将直接使用，否则根据条件判断
            mu_g: 气相粘度 (Pa·s)
            mu_l: 液相粘度 (Pa·s)
            temperature: 温度 (K)
            pressure: 压力 (Pa)
        Returns:
            float: 界面剪切力
        """
        # 计算相对速度
        v_r = v_g - v_l
        
        # 阶段2优化：添加速度饱和函数，避免二次增长
        v_ref = 1.0  # 特征速度参考值
        v_r_sat = np.tanh(abs(v_r)/v_ref) * np.sign(v_r)
        
        # 阶段2优化：添加指数阻尼因子，提高稳定性
        damping = np.exp(-0.5 * abs(v_r))
        
        # 改进的表面张力模型，考虑温度和压力影响
        sigma_0 = 0.0756 - 1.68e-4 * (temperature - 273.15)  # 基于温度的基础表面张力
        pressure_factor = 1.0 - 0.001 * (pressure / 1.0e6)  # 压力影响因子
        sigma = max(sigma_0 * pressure_factor, 0.001)  # 确保表面张力为正
        
        # 计算韦伯数
        We = rho_l * abs(v_r)**2 * d_h / max(sigma, 1e-10)
        
        # 界面面积密度(ai)和界面摩擦系数(fi)初始化
        a_i = 0.0  # 界面面积密度
        f_i = 0.0  # 界面摩擦系数
        
        # 如果未提供流型，使用流型识别器进行判断
        if flow_regime is None:
            # 阶段2优化：强制低空泡率使用气泡流，避免流型频繁切换
            if alpha < 0.1:
                flow_regime = 'bubbly'
            else:
                # 准备流型识别器所需的参数
                params = {
                    'diameter': d_h,
                    'rho_l': rho_l,
                    'rho_g': rho_g,
                    'mu_l': mu_l,
                    'mu_g': mu_g,
                    'sigma': sigma,  # 使用上面计算的表面张力
                    'pressure': pressure,
                    'temperature': temperature
                }
                
                # 调用流型识别器进行流型判断
                flow_regime = self.flow_regime_identifier.identify_regime(alpha, abs(v_r), **params)
        
        # 根据流型进行计算
        if flow_regime == 'bubbly':
            # 气泡流
            # 计算气泡直径，考虑表面张力和流动条件
            critical_diameter = 4.0 * sigma / (rho_l * abs(v_r)**2 + 0.1)
            bubble_diameter = min(0.05 * d_h, max(critical_diameter, 0.001))
            
            # 气泡雷诺数
            bubble_reynolds = abs(v_r) * bubble_diameter * rho_l / mu_l
            
            # 界面面积密度，考虑气泡流中气泡数量和分布
            a_i = 6.0 * alpha / max(bubble_diameter, 1e-6) * (1.0 + 0.5 * alpha)
            
            # 气泡拖曳系数计算，基于雷诺数和Eötvös数
            eotvos_number = abs(rho_l - rho_g) * 9.81 * bubble_diameter**2 / max(sigma, 1e-10)
            shape_factor = min(eotvos_number**0.5, 3.0)  # 气泡形状因子
            
            if bubble_reynolds < 1.0:
                # 斯托克斯区域
                C_d = 24.0 / max(bubble_reynolds, 0.1) * (1.0 + 0.1 * shape_factor)
            elif bubble_reynolds < 1000.0:
                # 过渡区域
                C_d = (24.0 / max(bubble_reynolds, 0.1)) * (1.0 + 0.15 * bubble_reynolds**0.687) * (1.0 + 0.1 * shape_factor)
            else:
                # 高雷诺数区域
                C_d = 0.44 * (1.0 + 0.15 * shape_factor)
            
            # 界面摩擦系数
            f_i = 0.75 * C_d * a_i
            
        elif flow_regime == 'slug':
            # 段塞流
            # 泰勒气泡特性长度(以水力直径为特征长度)
            taylor_bubble_length = d_h * (0.4 + 0.6 * alpha**0.5)
            
            # 界面面积密度，考虑单位体积内的界面面积
            a_i = 4.5 * alpha**0.3 * (1 - alpha)**0.5 / d_h
            
            # 界面摩擦系数，考虑段塞流的特性
            capillary_number = mu_l * abs(v_r) / max(sigma, 1e-10)  # 毛细管数(粘性力与表面张力比值)
            
            # 考虑表面张力对段塞流的影响
            st_correction = 1.0 + 0.8 * (sigma / max(rho_l * abs(v_r)**2 * taylor_bubble_length, 1e-10))**0.25
            
            # 计算段塞流界面摩擦系数
            f_i = 0.02 * (1.0 + 3.0 * alpha**0.5) * (1.0 + 50.0 * capillary_number**0.5) * st_correction * a_i
            
        elif flow_regime == 'churn':
            # 湍流流动
            # 界面面积密度，考虑高湍流度下的界面破碎增强
            a_i = 5.0 * alpha**0.5 * (1.0 - alpha)**0.3 / d_h
            
            # 计算湍流雷诺数
            turbulent_re = abs(v_r) * d_h * ((rho_g * alpha + rho_l * (1-alpha))) / ((mu_g * alpha + mu_l * (1-alpha)))
            
            # 湑流强度因子
            turbulent_intensity = min(0.2, 0.01 * (max(turbulent_re, 1.0) / 1000.0)**0.5)
            
            # 界面摩擦系数，考虑湍流增强和波动效应
            f_i = 0.035 * (1.0 + 10.0 * alpha * (1.0 - alpha)) * a_i * (1.0 + 5.0 * turbulent_intensity)
            
        else:  # annular
            # 环状流动
            # 液膜厚度计算
            film_thickness = d_h * (1.0 - alpha**0.5) / 4.0
            film_thickness = max(film_thickness, 0.0001 * d_h)  # 确保有最小液膜厚度
            
            # 界面面积密度，环状流中主要由中心气柱与环状液膜的界面决定
            a_i = 4.0 * alpha**0.25 / d_h
            
            # 计算无量纲参数
            film_reynolds = abs(v_l) * film_thickness * rho_l / mu_l  # 液膜雷诺数
            
            # 考虑界面波动的影响
            wave_factor = 1.0
            if film_reynolds > 1000:
                # 波动增强因子，液膜雷诺数大时波动增强
                wave_factor = 1.0 + 0.2 * (film_reynolds / 1000.0)**0.4
            
            # 环状流界面摩擦系数
            f_i = 0.005 * wave_factor * (1.0 + 300.0 * (film_thickness / d_h)) * a_i
        
        # 考虑子通道类型的影响
        f_i *= self.channel_coefficients[channel_type]['lateral_flow']
        
        # 优化：添加速度限制器，防止相对速度过大
        v_r_limited = np.clip(v_r, -5.0, 5.0)
        
        # 计算界面摩擦力（阶段2优化：使用饱和速度v_r_sat并乘以damping因子）
        raw_friction = f_i * rho_g * abs(v_r_sat) * v_r_sat * 0.05 * damping
        
        # 应用物理限制确保稳定性
        # 限制界面摩擦力不超过合理范围，防止出现数值不稳定
        max_drag = 0.1 * rho_l * abs(v_l)**2  # 基于液相动能的最大拖曳
        
        # 限制界面摩擦力的大小
        limited_friction = np.clip(raw_friction, -max_drag, max_drag)
        
        # 松弛因子应用，确保迭代稳定性（优化：从0.3降低到0.2）
        relaxation_factor = 0.2
        limited_friction *= relaxation_factor
        
        return limited_friction
        
    def calculate_lateral_interfacial_friction(self, alpha, rho_g, rho_l, v_g_lat, v_l_lat, d_h, channel_type=None, flow_regime=None, mu_g=1.8e-5, mu_l=1.0e-3, gap_width=None, temperature=298.15, pressure=1.0e5):
        """计算横向界面摩擦力
        
        Args:
            alpha: 空泡份额
            rho_g: 气相密度
            rho_l: 液相密度
            v_g_lat: 气相横向速度
            v_l_lat: 液相横向速度
            d_h: 水力直径
            channel_type: 子通道类型
            flow_regime: 流型字符串(可选)，如果提供将直接使用，否则根据条件判断
            mu_g: 气相粘度 (Pa·s)
            mu_l: 液相粘度 (Pa·s)
            gap_width: 间隙宽度 (m)
            temperature: 温度 (K)
            pressure: 压力 (Pa)
        
        Returns:
            float: 横向界面摩擦力
        """
        # 计算横向相对速度
        v_r_lat = v_g_lat - v_l_lat
        
        # 改进的表面张力模型，考虑温度和压力影响
        sigma_0 = 0.0756 - 1.68e-4 * (temperature - 273.15)  # 基于温度的基础表面张力
        pressure_factor = 1.0 - 0.001 * (pressure / 1.0e6)  # 压力影响因子
        sigma = max(sigma_0 * pressure_factor, 0.001)  # 确保表面张力为正
        
        # 子通道类型修正因子
        lateral_flow_factor = 1.0
        if channel_type is not None:
            lateral_flow_factor = self.channel_coefficients[channel_type]['lateral_flow']
        
        # 根据流型选择界面摩擦系数
        if flow_regime is None:
            # 准备流型识别器所需的参数
            params = {
                'diameter': d_h,
                'rho_l': rho_l,
                'rho_g': rho_g,
                'mu_l': mu_l,
                'mu_g': mu_g,
                'sigma': sigma,  # 使用上面计算的表面张力
                'pressure': pressure,
                'temperature': temperature
            }
            
            # 调用流型识别器进行流型判断（使用横向相对速度）
            flow_regime = self.flow_regime_identifier.identify_regime(alpha, abs(v_r_lat), **params)
        
        # 计算韦伯数
        We_lat = rho_l * abs(v_r_lat)**2 * d_h / max(sigma, 1e-10)
        
        # 计算界面摩擦系数（根据横向流动特性调整）
        C_i_base = self.C_i  # 基础界面摩擦系数
        
        # 根据流型调整界面摩擦系数
        if flow_regime == 'bubbly':
            # 气泡流 - 配合气泡变形和横向流动相互作用
            # 计算气泡直径
            bubble_diameter = min(0.03 * d_h, max(0.002, 4.0 * sigma / (rho_l * abs(v_r_lat)**2 + 1.0)))
            
            # 计算气泡雷诺数
            Re_bubble = rho_l * abs(v_r_lat) * bubble_diameter / mu_l
            
            # 气泡变形因子（横向流动时受表面张力的影响更大）
            deformation_factor = 1.0 + 0.15 * min(We_lat / 5.0, 1.0)
            
            # 横向流动的气泡流界面摩擦系数
            C_i = C_i_base * (1.0 + 20.0 * alpha * (1.0 - alpha)) * deformation_factor
            
        elif flow_regime == 'slug':
            # 弹状流 - 四周界面复杂性更高
            distortion_factor = 1.2 + 0.3 * We_lat**0.5
            C_i = C_i_base * distortion_factor * (1.0 + 2.0 * alpha)
            
        elif flow_regime == 'churn':
            # 搅拌流 - 使用更高的界面摩擦
            turbulence_factor = 1.5 * (1.0 + 0.5 * min(We_lat, 10.0) / 10.0)
            C_i = C_i_base * turbulence_factor * (1.0 + 5.0 * alpha * (1.0 - alpha))
            
        else:  # annular
            # 环状流 - 使用更低的界面摩擦
            film_factor = 0.8 * (1.0 + 0.2 * (1.0 - alpha)**0.5)
            C_i = C_i_base * film_factor
        
        # 应用子通道类型修正
        C_i *= lateral_flow_factor
        
        # 计算横向界面摩擦力
        lateral_interfacial_friction = C_i * rho_g * alpha * abs(v_r_lat) * v_r_lat
        
        # 同过数值阻尼提高计算稳定性
        damping_factor = np.exp(-0.5 * abs(v_r_lat))  # 速度大时慢慢减小阻尼
        damped_friction = lateral_interfacial_friction * damping_factor
        
        # 限制摩擦力范围，避免非物理数值
        max_friction = 0.2 * rho_l * abs(v_l_lat)**2  # 基于液相动能限制
        limited_friction = np.clip(damped_friction, -max_friction, max_friction)
        
        return limited_friction
        
    def calculate_turbulent_mixing_momentum(self, alpha, rho_g, rho_l, v_g, v_l, phase, d_h, channel_type=None, flow_regime=None, mu_g=1.8e-5, mu_l=1.0e-3, temperature=298.15, pressure=1.0e5, mixing_area=None):
        """计算湍流混合动量
        
        基于流型的湍流混合模型，考虑各种流型特性对湍流混合的影响，
        包括气泡流、段塞流、搅拌流和环状流等情况下的不同湍流混合机制。
        
        Args:
            alpha: 空泡份额
            rho_g: 气相密度
            rho_l: 液相密度
            v_g: 气相轴向速度
            v_l: 液相轴向速度
            phase: 相态（'vapor'或'liquid'）
            d_h: 水力直径
            channel_type: 子通道类型
            flow_regime: 流型字符串(可选)，如果提供将直接使用，否则根据条件判断
            mu_g: 气相粘度 (Pa·s)
            mu_l: 液相粘度 (Pa·s)
            temperature: 温度 (K)
            pressure: 压力 (Pa)
            mixing_area: 混合区域面积 (m²)，可选，若未提供则估算
        
        Returns:
            float: 湍流混合动量
        """
        # 计算相对速度，这是湍流混合强度的主要驱动因素
        v_r = v_g - v_l
        
        # 根据相态确定物性参数
        if phase == 'vapor':
            rho = rho_g
            mu = mu_g
            phase_fraction = alpha
            vel = v_g
        else:  # 'liquid'
            rho = rho_l
            mu = mu_l
            phase_fraction = 1.0 - alpha
            vel = v_l
            
        # 限制相分数在有效范围内，避免数值问题
        phase_fraction = np.clip(phase_fraction, 1e-6, 1.0-1e-6)
        
        # 改进的表面张力模型，考虑温度和压力影响
        sigma_0 = 0.0756 - 1.68e-4 * (temperature - 273.15)  # 基于温度的基础表面张力
        pressure_factor = 1.0 - 0.001 * (pressure / 1.0e6)    # 压力影响因子
        sigma = max(sigma_0 * pressure_factor, 0.001)         # 确保表面张力为正
        
        # 子通道类型修正因子
        turbulent_mixing_factor = 1.0
        if channel_type is not None:
            turbulent_mixing_factor = self.channel_coefficients[channel_type]['turbulent_mixing']
        
        # 如果未提供流型，使用流型识别器进行判断
        if flow_regime is None:
            # 准备流型识别器所需的参数
            params = {
                'diameter': d_h,
                'rho_l': rho_l,
                'rho_g': rho_g,
                'mu_l': mu_l,
                'mu_g': mu_g,
                'sigma': sigma,  # 使用上面计算的表面张力
                'pressure': pressure,
                'temperature': temperature
            }
            
            # 调用流型识别器进行流型判断
            flow_regime = self.flow_regime_identifier.identify_regime(alpha, abs(v_r), **params)
        
        # 湍流混合系数基础值
        beta_base = self.BETA  # 基础湍流交混系数
        
        # 根据流型计算湍流混合系数
        if flow_regime == 'bubbly':
            # 气泡流 - 气泡诱导湍流增强
            # 计算气泡直径，考虑表面张力和流动条件
            bubble_diameter = min(0.03 * d_h, max(0.002, 4.0 * sigma / (rho_l * abs(v_r)**2 + 1.0)))
            
            # 气泡雷诺数
            bubble_reynolds = max(abs(v_r) * bubble_diameter * rho_l / mu_l, 0.1)
            
            # Eötvös数(浮力与表面张力的比值)
            eotvos_number = abs(rho_l - rho_g) * 9.81 * bubble_diameter**2 / sigma
            
            # 气泡诱导湍流因子
            bubble_induced_factor = 1.0 + 0.8 * min(eotvos_number**0.25, 3.0) * (alpha**0.5)
            
            # 气泡湍流强度因子 - 随雷诺数增加而增强
            if bubble_reynolds < 100:
                # 低雷诺数区域 - 层流主导
                turbulence_intensity = 0.02 * (bubble_reynolds / 100)**0.5
            elif bubble_reynolds < 1000:
                # 过渡区域
                turbulence_intensity = 0.02 + 0.08 * ((bubble_reynolds - 100) / 900)**0.7
            else:
                # 高雷诺数区域 - 完全湍流
                turbulence_intensity = 0.1 + 0.05 * min((bubble_reynolds / 1000)**0.2, 2.0)
            
            # 气泡流湍流混合系数
            beta = beta_base * bubble_induced_factor * (1.0 + 5.0 * turbulence_intensity) * (1.0 + 10.0 * alpha * (1.0 - alpha))
            
        elif flow_regime == 'slug':
            # 段塞流 - 泰勒气泡和液体段落交替引起的湍流
            # 泰勒气泡特性长度(以水力直径为特征长度)
            taylor_bubble_length = d_h * (0.4 + 0.6 * alpha**0.5)
            
            # 段塞流雷诺数 - 使用混合密度
            slug_reynolds = abs(v_r) * d_h * (rho_g * alpha + rho_l * (1-alpha)) / (mu_g * alpha + mu_l * (1-alpha))
            
            # 计算界面波动影响因子 - 段塞流中相界面波动明显
            wave_factor = 1.0 + 0.5 * min((slug_reynolds / 1000)**0.4, 2.0)
            
            # 段塞流湍流混合系数
            beta = beta_base * wave_factor * (1.0 + 3.0 * alpha**0.4 * (1.0 - alpha)**0.4)
            
        elif flow_regime == 'churn':
            # 湍流流动 - 高度湍流混合区域
            # 计算湍流雷诺数 - 使用混合物性
            turbulent_re = abs(v_r) * d_h * (rho_g * alpha + rho_l * (1-alpha)) / (mu_g * alpha + mu_l * (1-alpha))
            
            # 湍流强度因子 - 搅拌流中湍流强度更高
            turbulent_intensity = min(0.2, 0.02 * (turbulent_re / 1000.0)**0.6)
            
            # 湍流混合系数 - 搅拌流中湍流混合最强
            beta = beta_base * (2.0 + 5.0 * alpha * (1.0 - alpha)) * (1.0 + 7.0 * turbulent_intensity)
            
        else:  # annular
            # 环状流 - 液膜和气芯之间的湍流混合
            # 液膜厚度计算
            film_thickness = d_h * (1.0 - alpha**0.5) / 4.0
            film_thickness = max(film_thickness, 0.0001 * d_h)  # 确保最小液膜厚度
            
            # 液膜雷诺数
            film_reynolds = abs(v_l) * (2 * film_thickness) * rho_l / mu_l
            
            # 界面波动因子 - 环状流中界面波动较弱
            wave_factor = 1.0
            if film_reynolds > 1000:
                wave_factor = 1.0 + 0.1 * min((film_reynolds / 1000.0)**0.3, 2.0)
            
            # 环状流湍流混合系数 - 随液膜雷诺数增加而增强
            beta = beta_base * wave_factor * (0.6 + 0.4 * (1.0 - alpha)**0.5)
        
        # 应用子通道类型修正
        beta *= turbulent_mixing_factor
        
        # 相对速度效应 - 相对速度增大时湍流混合增强
        beta *= (1.0 + 0.2 * abs(v_r) / max(abs(vel), 0.01))
        
        # 估算湍流混合面积（如果未提供）
        if mixing_area is None:
            # 简单估计为水力直径周长与长度的乘积的一部分
            mixing_area = np.pi * d_h * d_h * 0.25  # 假设为圆形横截面积
        
        # 计算湍流混合动量
        mixing_momentum = 0.0
        if abs(v_r) > 1e-5:  # 避免除数过小引起数值问题
            # 使用模拟湍流黏度的方法
            artificial_eddy_viscosity = beta * d_h * abs(vel)
            mixing_momentum = artificial_eddy_viscosity * rho * v_r * mixing_area
            
            # 数值阻尼处理，提高计算稳定性
            damping_factor = np.exp(-0.3 * abs(v_r))  # 速度差大时提供阻尼
            mixing_momentum *= damping_factor
        
        # 物理上限制混合动量大小，避免非物理结果
        max_allowed_momentum = 0.3 * rho * vel**2 * mixing_area
        mixing_momentum = np.clip(mixing_momentum, -max_allowed_momentum, max_allowed_momentum)
        
        return mixing_momentum

    def calculate_lateral_momentum_exchange(self, alpha, rho_g, rho_l, v_g, v_l, v_g_adj, v_l_adj, gap_length, gap_width):
        """计算横向动量交换
        
        Args:
            alpha: 空泡份额
            rho_g: 气相密度
            rho_l: 液相密度
            v_g: 本通道气相速度
            v_l: 本通道液相速度
            v_g_adj: 相邻通道气相速度
            v_l_adj: 相邻通道液相速度
            gap_length: 间隙长度
            gap_width: 间隙宽度
        """
        # 计算横向速度差
        dv_g = v_g - v_g_adj
        dv_l = v_l - v_l_adj
        
        # 计算横向流动面积
        A_lat = gap_length * gap_width
        
        # 计算横向动量交换（考虑FTM）
        M_lat_g = -self.FTM * rho_g * alpha * A_lat * dv_g * abs(dv_g)
        M_lat_l = -self.FTM * rho_l * (1-alpha) * A_lat * dv_l * abs(dv_l)
        
        return M_lat_g, M_lat_l
        
